Spatial Analysis

Methods

Spatial data preparation

Load data

Maps data

We use three different shapefiles for the continental U.S. land mass, the State waters of maine, new hampshire, massachusets, connecticut, rhode island, new york, new jersey, delaware, maryland, virginia, north carolina and the North East EEZ. Moreover, we use the NOAA dataset of landings per port.

1.- Land shapefile

Covers the US land territory for visualization. Data provided from the map package.

States <- c("maine", "new hampshire", "massachusetts", "connecticut", "rhode island", "new york", "new jersey", "delaware", "maryland", "virginia", "north carolina") 

# -------------- #
# US State Map (land)
# -------------- #

land_sf <- st_as_sf(map("state", plot = FALSE, fill = TRUE)) %>% 
  filter(ID %in% States) %>% 
  mutate(abrev = c("CT","DE","ME","MD","MA","NH","NJ","NY","NC","RI","VA"),
         state = str_to_sentence(ID)
         )

# Visual exploration (o.k! )
ggplot(land_sf) +
  geom_sf(aes(fill = state)) +
  scale_fill_manual(values = state_pallet) +
  theme_dark()

NA
NA

2.- EEZ shapefile

Used the Sea Around Us shapefle updated to June 2016.

# US EEZ map
path_world <- my_path("G", extra_path = "Spatial/SAU/SAU_Shapefile", name = "SAUEEZ_July2015.shp")

# Load it!
eez_sf <- st_read(path_world) %>%
  rename(eez_name = Name) %>%
  st_transform(crs = 4326) %>% # 4326
  filter(eez_name == "USA (East Coast)") %>%
  st_simplify(preserveTopology = TRUE, dTolerance = 0.1) #0.1 for paper

Reading layer `SAUEEZ_July2015' from data source 
  `/Volumes/Enterprise/Data/Spatial/SAU/SAU_Shapefile/SAUEEZ_July2015.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 280 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -63.66443 xmax: 180 ymax: 87.02394
Geodetic CRS:  WGS 84

Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance) :
  st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

# Visual exploration (o.k! )
ggplot(eez_sf) +
  geom_sf() +
  theme_dark()

3.- State waters

Covers the state waters of the NE US states. Data from data.gov.

License: No license information was provided. If this work was prepared by an officer or employee of the United States government as part of that person’s official duties it is considered a U.S. Government Work.

state_sf <-  st_read("~/OneDrive - UW-Madison/FederalAndStateWaters/FederalAndStateWaters.gdb") %>%
  janitor::clean_names() %>% 
  mutate(jurisdicti = str_to_lower(jurisdicti)) %>% 
  filter(jurisdicti %in% States) %>% 
  mutate(state = str_to_sentence(jurisdicti)) %>% 
  st_transform(crs = 4326) %>%   # for matching projections
  st_simplify(preserveTopology = TRUE, dTolerance = 10000) #0.1 for paper

Reading layer `FederalAndStateWaters' from data source 
  `/Users/jepa/OneDrive - UW-Madison/FederalAndStateWaters/FederalAndStateWaters.gdb' 
  using driver `OpenFileGDB'

Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet,  :
  GDAL Message 1: organizePolygons() received a polygon with more than 100 parts. The processing may be really slow.  You can skip the processing by setting METHOD=SKIP, or only make it analyze counter-clock wise parts by setting METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock wise defined

Simple feature collection with 40 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -20037510 ymin: -1972646 xmax: 20037510 ymax: 12766910
Projected CRS: WGS 84 / World Mercator

Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance) :
  st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

# Visual exploration (o.k! )
ggplot(state_sf) +
  geom_sf(aes(fill = state)) +
  scale_fill_manual(values = state_pallet) +
  theme_dark()

4.- Fishing ports

4.1 NOAA’s data

Dataset provider: Joe Myers Joseph.Myers@accsp.org NOAA

Non-confidential dataset containing annual landings/value summaries by port for black sea bass, scup and summer flounder, from NC-ME for 1970-2020.

In order to make the report non-confidential, I am displaying as much of the port-level summary data as I can – and for cases where the port-level was confidential, I have first attempted to roll-up any confidential port-level summaries into state-level summaries within the same year. If those state-level summaries continue to be confidential, I have further rolled-up those confidential ports across states into a coast-level summary. Even at the coast level, there were several cases where not enough confidential ports could be combined to produce a non-confidential summary – and so for those rows it was necessary to redact the final landings/value data. Fortunately, this only amounts to less than 0.5% of the total landings data (characterized in the “Pivot by Confidentiality” tab of the spreadsheet).

License???.

4.1.1 Top ports

# Get number of ports in the NOAA database
overall_port_data <- ports_data %>%
  filter(year >= 2010 & year <= 2019) %>% 
  group_by(species_name,port_name,state_postal) %>% 
  summarise(
    year_ton = sum(total_live_pounds),
    year_value = sum(total_value)
  ) 

`summarise()` has grouped output by 'species_name', 'port_name'. You can override using the `.groups` argument.

4.1.2 Ports coordinates

4.1.1 Visualization

Create base shapefile for computation

As a first step we need to divide the NE US EEZ among the different states. For that, we expanded a buffer-polygon to the 200 nautical miles to then estimate the percentage that each expanded polygon occupied. For the initial buffer point we used two approaches; state waters and US ports. Note that in all cases these areas overlapped and percentages accounted for it. We did this by following these steps:

• 1. Expand state waters / US ports using a buffer
• 2. Grid that buffer on a 0.3 resolution
• 3. Crop the buffer to the EEZ

1. Expand Spatial regions (a.k.a buffers)

1.1 State Waters

We set a buffer of 410000 m (410 km, ~ 221 nm) that overshoots the EEZ a bit, but is eventually cropped

# Buffer state waters
state_bf = st_buffer(state_sf, 4) %>% 
  st_transform(4326) %>% # to match shape
  st_set_crs(4326)

# ------------------------------ #
# Testing and visualizing the buffer 
# ------------------------------ #

state_bf_plot <- ggplot() +
  geom_sf(data = state_sf, aes(color = state), fill = NA) +
  geom_sf(data = eez_sf, aes(), fill = NA) +
  geom_sf(data = state_bf, aes(color = state), fill = NA) +
  scale_color_manual(values = state_pallet)+
  theme_dark()
#   

# ggsave("../Results/Partial/state_waters_buffer.png",state_bf_plot)

# ------------------------------ #

state_bf_plot

1.1 US. Ports

We set a buffer of 5 deg (510 km, ~ 221 nm) that overshoots the EEZ a bit, but is eventually cropped

2. Expand grid within buffer

Here we expand a grid within the buffer so we can estimate the proportion of each state.

Note: Dark grey shaded area is the grid and is the same for both approaches

2.1 Create grid for interpolation

2.1 Merge grid and buffers

Once we have a gridded area, we converted the grid to a sf so we can merge it with the buffered states and ports and finally filter out everything outside the states polygon

2.1.1 State Waters

# ---------------- #
# Convert grid to sf
# ---------------- #
grid_sw_sf <- st_as_sf(ne_grid,
             coords = c("lon", "lat"),
             crs = 4326) %>% 
  st_join(state_bf) %>% 
  filter(!is.na(state))

# Create data frame for future computations
# Note, will be used in next chunk
grid_sw_bf_dt <- as.data.frame(grid_sw_sf) %>%
    select(index,state)
    # group_by(state) %>% 
    # summarise(length(index))


# ---------------- #
# [TEST] 
# Visualization of grid
# ---------------- #

gridExtra::grid.arrange(
  # Overall (overlapping) position
  ggplot(grid_sw_sf) +
    geom_sf(data = state_sf, aes(), fill = NA) +
    geom_sf(aes(color = state), alpha = 0.3) +
    geom_sf(data = eez_sf, aes(),fill =NA) + 
    scale_color_manual(values = state_pallet) +
    theme_dark() +
    theme(legend.position = ""),
  # Showing each state separatley
  ggplot() +
    geom_sf(data = grid_sw_sf, aes(color = state),size = 0.1, alpha = 0.5) +
    facet_wrap(~state) +
    theme_dark() +
    theme(legend.position = "top") +
    scale_color_manual(values = state_pallet),
  nrow = 1) 

2.1.1 Fishing Ports

3. Crop buffers to EEZ

Finally, we crop the grided buffers to within the EEZ to capture the actual water space.

Note: This step takes quite a while because of the size of the EEZ shapefile. No, you cannot use st_simplify() here

3.1 State Waters

grid_eez_sw_sf <- st_intersection(grid_sw_sf,eez_sf)

Error in st_intersection(grid_sw_sf, eez_sf) : 
  object 'grid_sw_sf' not found

3.2 Fishing Ports

3.3 Differences at croped level

By looking at the plots we do not see that the number of grid cells that each state has is slightly different depending on the approach. This comes to light when we count the number of gridcells that each state gets in each approach. We see that the state waters one is substantially higher.

We can visualize the difference of the cropped area for each state next. Again, the fishing port approach results in less area for all states, but such area is not proportional among states. For example, Delaware looses roughly 35% of fishing grounds with this approach

grid_eez_sw_df %>% 
  mutate(approach = "state_waters") %>% 
  bind_rows(grid_eez_fp_df) %>%
  mutate(approach = ifelse(is.na(approach),"fishing_port",approach)) %>% 
  group_by(state,approach) %>% 
  summarise(n_grid = n()) %>% 
  spread(approach,n_grid) %>% 
  mutate(difference =state_waters-fishing_port,
         percentage_sw = round((difference/state_waters)*100)
         )

And here is the map visualization of the cropped difference

grid_eez_sw_df %>% 
  mutate(approach = "state_waters") %>% 
  bind_rows(grid_eez_fp_df) %>%
  mutate(approach = ifelse(is.na(approach),"fishing_port",approach)) %>% 
  ggplot() +
  geom_tile(
    aes(
      x = lon,
      y = lat,
      fill = approach
    ),
    alpha = 0.5
  ) +
  facet_wrap(~state) +
  scale_color_brewer(palette = "Set2") +
  scale_fill_brewer(palette = "Set2") +
  theme_dark()

Interpolation rutine

In this step we interpolate the survey data within the previously created grid following a Triangular Irregular Surface method.

Functions needed

We need to create a couple of functions to run the whole process

Interpolation main fx (tis).

This is the main function used to interpolate the data per year. It follows a Triangular Irregular Surface method using the interp::interp() function. If you want to see the function clic on code

tis <- function(input_data, grid, yr, taxa, reg){
  
  # --------------- #
  # Testing
  # print(paste(yr))
  # yr = 1976
  # --------------- #
  
  # Filter data
  data <- input_data %>% 
    filter(year == yr,
           spp == taxa,
           region == reg
    ) %>% 
    # Average duplicated hauls in the same spot
    group_by(region,year,spp,lat,lon) %>%
    summarise_at(vars(wtcpue),
                 mean,na.rm = T)
    
  # Only interpolate cases where there is more than 3 rows
  # Marked by the function 
    if(nrow(data) <= 3){
      fit_tin <- tibble()
    }else{
      
      # Triangular Irregular Surface
      fit_tin <- interp::interp( # using {interp}
        x = data$lon,           # the function actually accepts coordinate vectors
        y = data$lat,
        z = data$wtcpue,
        xo = grid$lon,     # here we already define the target grid
        yo = grid$lat,
        output = "points"
      ) %>% 
        bind_cols() %>% 
        bind_cols(grid) %>%
        mutate(year = yr,
               region = reg,
               spp = taxa) %>% 
        select(index, 
               # state,
               year, region, spp, lon=x, lat=y, value = z)
      
    }
   
    return(fit_tin)
}

# Test me
# Needs variables in Control panel
# Test no data: "Alosa aestivalis", reg = "Northeast US Fall", yr = 1974 
# tis(input_data = ocean_data, grid = grid_eez_df, taxa = "Illex illecebrosus", reg = "Northeast US Fall", yr = 1973)



# lapply(years,tis,input_data = ocean_data, grid = grid_eez_df, taxa = "Illex illecebrosus", reg = regions[2])

Run function

This is a sub-function that runs the tis() function by taxa and region. It saves the output as a .csv file for each species.

run_tis <- function(input_data, grid, years, taxa, reg){
  
  
  # Run tis for species and surveys
  for(r in 1:2){
    
    partial_df <- bind_rows(
      lapply(years,tis,input_data = input_data, grid = grid, taxa = taxa, reg = reg[r])
    )
    
    if(r == 1){
      historic_tif <- partial_df
    }else{
      historic_tif <- bind_rows(historic_tif,partial_df)
    }
    
  }
  
  # ----------------------- #
  # Save dataset per species
  # ----------------------- #
  
  # Set file name
  name <- paste0("tif_",str_replace(taxa," ","_"),".csv")
  
  complement <- paste0("Partial/Interpolation/")
  
  # Set path name
  save_path <- my_path("R",complement)
  
  # Set complete path
  save_name <- paste0(save_path,name)
  
  # Create folder if it does not exist
  if(file.exists(save_path) == F){
    dir.create(save_path)
  }
  
  #  Save file
  write_csv(historic_tif,save_name)
  
  return_msg <- paste("Interpolation done for", taxa)
  
  return(return_msg)
  
  
}

# Test me
    # run_tis(input_data = ocean_data, grid = grid_eez_fp_df, taxa = "Centropristis striata", years = seq(1970,1980), reg = regions, method = "testa")

Survey data

The interpolation was done with NOAA Northeast Fisheries Science Center Spring and Fall Bottom Trawl Surveys data provided by Ocean adapt. Data was accessed trough the Github.

In primary publications using data from the database, please cite Pinsky et al. 2013. Marine taxa track local climate velocities. Science 341: 1239-1242 doi: 10.1126/science.1239352, as well as the original data sources.

• Using only the Northeast US Fall and Spring bottom trawl survey data for now

Splitting up data

• No part on spatial analysis. Can be ignored.

This is just a sub-step to split up the data into single species files. This makes the app faster as it only needs to load species specific data, rather than all the data at de beginning.

Control Pannel

This is where we load the required data and prepare to run the interpolation function. Note that some of the data here has been previously created

# Load grid df
# For fishing ports
inter_grid <- my_path("D","Spatial","inter_grid_df.csv", read = T)

# Run interpolation for all years
years <- seq(1973,2019,1)

# regions
regions <- c("Northeast US Fall" , "Northeast US Spring")

# species list
spp <- ocean_data %>% 
  filter(region %in% regions,
         spp != "NA") %>% 
  pull(spp) %>% 
  unique()


# Double check runs

spp_runs <- tibble(taxa = (list.files(my_path("R","Partial/Interpolation/fp")))) %>%
  mutate(
    taxa = str_remove(taxa, "tif_"),
    taxa = str_remove(taxa, ".csv"),
    taxa = str_replace(taxa, "_", " ")
  ) 

# Missing runs
spp <- tibble(taxa=spp) %>% 
  anti_join(spp_runs) %>% 
  pull(taxa)

Run routine

Here we just run the routine for each of the species present in the Northeast US Fall and Spring surveys between 1973 and 2019.

• Note there are 43 identified taxa
• Some taxa do not have presence data in some years

# # single species run
run_tis(input_data = ocean_data,
        grid = grid_eez_fp_df,
        years = years,
        reg = regions,
        taxa = "Stenotomus chrysops",
        method = "fp"
)


# Run them all in parallel
lapply(Species,
       run_tis, 
       input_data = ocean_data,
       grid = inter_grid,
       years = years,
       reg = regions
)

Results

Results are now for eight species managed under Mid-Atlantic Council Management Plans according to NOAA.

State allocations of the commercial black sea bass coastwide quota were originally implemented in 2003 as part of Amendment 13, loosely based on historical landings from 1980- 2001.

Scomber scombrus, Peprilus triacanthus (butterfish), Illex illecebrosus (shortfin squid), Paralichthys dentatus (summer flounder), Stenotomus chrysops (Scoop), Centropristis striata (black sea bass), Pomatomus saltatrix (Bluefish), Lopholatilus chamaeleonticeps (Golden tilefish), Caulolatilus microps (blueline tilefish) and Clupea harengus

# Spatial data
States <- c("maine", "new hampshire", "massachusetts", "connecticut", "rhode island", "new york", "new jersey", "delaware", "maryland", "virginia", "north carolina") 

# US State Map (land)
land_sf <- st_as_sf(map("state", plot = FALSE, fill = TRUE)) %>% 
  filter(ID %in% States) %>% 
  mutate(abrev = c("CT","DE","ME","MD","MA","NH","NJ","NY","NC","RI","VA"),
         state = str_to_sentence(ID)
         )

# US EEZ map
path_world <- my_path("G", extra_path = "Spatial/SAU/SAU_Shapefile", name = "SAUEEZ_July2015.shp")


### Previouslly created grids 

grids <- my_path("D", "Spatial/grid_eez_fp", name = "grid_eez_fp_df.csv", read = T) %>%
  mutate(
    spatial = "fp"
  ) %>% 
  bind_rows(
    my_path("D", "Spatial/grid_eez_sw", name = "grid_eez_sw_df.csv", read = T)
  ) %>% 
  select(-spp) %>% 
  mutate(
    spatial = ifelse(is.na(spatial),"sw",spatial)
  ) %>% 
  distinct(.keep_all = T)
  

# Species data

# Managed species by Mid-Atlantic Council Management Plans, and State waters
# https://www.fisheries.noaa.gov/new-england-mid-atlantic/commercial-fishing/new-england-mid-atlantic-fishery-management-plans

Species <- c(
  "Paralichthys dentatus", #summer flounder
  "Stenotomus chrysops", #Scup"
  "Centropristis striata"#, # black sea bass
  )

spp_paths <- list.files(my_path("R",extra_path = "Partial/Interpolation/"))

# Some gibiris for names
taxon_list <- gsub("_"," ",spp_paths)
taxon_list <- gsub("tif ","",taxon_list)
taxon_list <- gsub(".csv","",taxon_list)

taxon_list <- taxon_list[taxon_list %in% Species]

# Now set a list of files to read
taxon_read <- paste0("tif_",taxon_list,".csv")
taxon_read <- gsub(" ","_",taxon_read)

taxon_read <- spp_paths[spp_paths %in% taxon_read]

#### Load sw interpolation 

taxon_sw <- my_path("R",extra_path = "Partial/Interpolation/", name = taxon_read)

# Load all sp in one table
historic_spp <- bind_rows(lapply(taxon_sw, fread)) %>% 
  left_join(grids,
            by = c("index","lon","lat")
  ) %>%
  filter(!is.na(spatial)) %>% 
  mutate(
    season = ifelse(str_detect(region,"Spring"),"Spring","Fall"),
    label = ifelse(spp == "Paralichthys dentatus" & year >= 1980 & year <= 1986,"Reference",
                   ifelse(spp == "Stenotomus chrysops" & year >= 1988 & year <= 1992,"Reference",
                          ifelse( spp =="Centropristis striata" & year >= 1980 & year <= 2001,"Reference",
                                 ifelse(year > 2001,"Today",NA)
                          )
                   )
    )
  )

# historic_spp %>% 
#   filter(is.na(spatial)) %>% 
#   filter(year == 1973,
#          region == "Northeast US Fall") %>% 
#   ggplot() +
#     geom_tile(
#       aes(
#         x = lon,
#         y = lat,
#         fill = index
#       )
#     ) +
#   facet_grid(spp~state) + 
#   geom_sf(data = land_sf,aes()) +
#   geom_sf(data = eez_sf, aes(), fill =NA)

unique(historic_spp$spp)
unique(historic_spp$region)
unique(historic_spp$season)

# Periods
periods <-tibble(
  order = c(rep("a",12),
            rep("b",12),
            rep("c",12),
            rep("d",11)
  ),
  label = c(rep("1973-1984",12),
            rep("1985-1996",12),
            rep("1997-2008",12),
            rep("2009-2019",11)
  ),
  year = c(seq(1973,1984,1),
           seq(1985,1996,1),
           seq(1997,2008,1),
           seq(2009,2019,1)
  )
  )


state_lat <- historic_spp %>%
  group_by(state) %>% 
  summarise(
    order = mean(lat)
  ) %>% 
  filter(!is.na(state)) %>% 
  mutate(abrev = c("CT","DE","ME","MD","MA","NH","NJ","NY","NC","RI","VA")
  ) %>% 
  arrange(desc(order)) %>% 
  mutate(lat = order,
    # order = letters[1:11],
    order = c("b","a",letters[3:11]),)



# State pallet

state_pallet <- c(wes_palette(n = 5, name = "Darjeeling1"),
                               wes_palette(n = 5, name = "Darjeeling2"),
                               # wes_palette(n = 2, name = "GrandBudapest1")
                  "#FD6467"
                  )



# print for notebook
head(historic_spp)

Regulatory units result

#### Data needed ####

# Get state order from North to south
state_order <- my_path("D","Spatial/grid_eez_sw", name = "grid_eez_sw_df.csv",read = T) %>%
  group_by(state) %>% 
  summarise(
    order = mean(lat)
  ) %>% 
  mutate(abrev = c("CT","DE","ME","MD","MA","NH","NJ","NY","NC","RI","VA")
  ) %>% 
  arrange(desc(order)) %>% 
  mutate(order= c("b","a",letters[3:11])) %>% # weirdly maine goes below
  arrange(order)



# Load grid of state waters
grid_eez_sw_sf <-  st_read(my_path("D","Spatial/grid_eez_sw", name = "grid_eez_sw_sf.shp")) %>%
  select(index,state,geometry) %>% 
  mutate(method = "state_waters",
         state = str_to_sentence(state)) %>% 
  left_join(state_order)

grid_eez_sw_sf$group <- factor(grid_eez_sw_sf$abrev,      # Reordering group factor levels
                               levels = state_order$abrev)


# Load grid of fishing ports
grid_eez_fp_sf <-  st_read(my_path("D","Spatial/grid_eez_fp", name = "grid_eez_fp.shp")) %>%
  select(index,abrev=lndng_s,geometry) %>% 
  mutate(method = "fishing_port"#,
         # state = str_to_sentence(landing_port)
         )%>% 
  left_join(state_order)


grid_eez_fp_sf$group <- factor(grid_eez_fp_sf$abrev,      # Reordering group factor levels
                               levels = state_order$abrev)

grid_eez_fp_sf <- arrange(grid_eez_fp_sf,group)
# Plotting 
# Make sure st_simplify() is run for eez sf

# State waters plot

p_sw <- gridExtra::grid.arrange(
  # Overall (overlapping) position
  ggplot(grid_eez_sw_sf) +
    geom_sf(data = eez_sf, aes(), fill = "white") + 
    geom_sf(aes(color =group , fill = group), alpha = 0.3) +
    geom_sf(data = land_sf, aes()) +
    ggsflabel::geom_sf_label_repel(data = land_sf, aes(label= abrev), 
                                   size = 4,
                                   box.padding = 0.10,
                                   hjust = 1) +
    scale_color_manual(values = state_pallet) +
    scale_fill_manual(values = state_pallet) +
    my_ggtheme_p(leg_pos = "",
                 ax_tx_s = 13) +
    coord_sf(ylim = c(30,48)) +
    scale_y_continuous(breaks = c(30,35,40,45))+
    labs(x = "", y = "", title = "A) State waters approach") +
    theme(plot.title = element_text(size = 20)),
  # Showing each state separately
  ggplot(grid_eez_sw_sf) +
    geom_sf(data = land_sf, aes(), fill = "grey80") +
    geom_sf(aes(color = group),size = 0.1, alpha = 0.3) +
    facet_wrap(~ group) +
    theme(legend.position = "top") +
    scale_color_manual(values = state_pallet,
                       # labels = unique(grid_eez_sw_sf$group)) +
                       labels = grid_eez_sw_sf %>% arrange(order) %>%  pull(state) %>% unique()) +
    ggtitle("") +
    my_ggtheme_p(leg_pos = "",
                 ax_tx_s = 11,
                 axx_tx_ang = 45,
                 hjust = 1
    ),
  nrow = 1)


##### ------------- ####
# Fishing Ports plot
##### ------------- ####


grid_eez_fp_sf <-
  grid_eez_sw_sf %>% 
  filter(!group %in% grid_eez_fp_sf$group) %>% 
  mutate(index = NA) %>% 
  bind_rows(grid_eez_fp_sf)
  
  


p_fp <- gridExtra::grid.arrange(
  # Overall (overlapping) position
  ggplot(data = subset(grid_eez_fp_sf, index != "NA")) +
    geom_sf(data = eez_sf, aes(), fill = "white") + 
    geom_sf(aes(color = group, fill = abrev), alpha = 0.3) +
    geom_sf(data = land_sf, aes()) +
    ggsflabel::geom_sf_label_repel(data = land_sf, aes(label= abrev), 
                                   size = 4,
                                   box.padding = 0.10,
                                   hjust = 1) +
    # scale_color_manual(values = state_pallet) +
    # scale_fill_manual(values = state_pallet) +
    scale_color_manual(values = c(state_pallet[3:4],state_pallet[6:13])) +
    scale_fill_manual(values = c(state_pallet[3:4],state_pallet[6:13])) +
    my_ggtheme_p(leg_pos = "",
                 ax_tx_s = 13) +
    coord_sf(ylim = c(30,48)) +
    scale_y_continuous(breaks = c(30,35,40,45)) +
    labs(x = "", y = "", title = "B) Fishing ports approach") +
    theme(plot.title = element_text(size = 20)),
  # Showing each state separately
  ggplot() +
    geom_sf(data = subset(grid_eez_fp_sf, index = "NA"), color = "white") +
    geom_sf(data = subset(grid_eez_fp_sf, index != "NA"), aes(color = group),size = 0.1, alpha = 0.3) +
    geom_sf(data = land_sf, aes(), fill = "grey80") +
    facet_wrap(~ group) +
    geom_point(data = rename(ports_df, group = landing_state), aes(x = lon, y = lat), color = "#E6A0C4") +
    theme(legend.position = "top") +
    # scale_color_manual(values = state_pallet,
    scale_color_manual(values = c(state_pallet[3:4],state_pallet[6:13]),
                       # labels = unique(grid_eez_fp_sf$group)) +
                       labels = grid_eez_fp_sf %>% arrange(order) %>%  pull(state) %>% unique()) +
    ggtitle("") +
    my_ggtheme_p(leg_pos = "",
                 ax_tx_s = 11,
                 axx_tx_ang = 45,
                 hjust = 1
    ),
  nrow = 1)



combined_plot <- gridExtra::grid.arrange(p_sw,p_fp,
                                         bottom = "Longitude",
                                         left = "Latitude")

ggsave(filename = "buffer_figure.jpg",
       plot = combined_plot,
       path = my_path("R","Figures"),
       width = 11,
       height = 11
       )

Interpolation result

# Interpolation data
interpol_result <- historic_spp %>% 
  filter(!is.na(label)) %>% 
  group_by(label,index,spp,season,lat,lon) %>% 
  summarise(mean_cpue = mean(value,na.rm=T)) %>% 
  filter(!is.na(mean_cpue))


# Paper numbers

# interpol_result %>% 
#   filter(mean_cpue >0) %>% 
#   View()

# --------------------- #
# Biomass Centroid shift
# --------------------- #

# Estimate total interpolation
total_inter <- historic_spp %>% 
  filter(!is.na(label)) %>% 
  group_by(spp,index,season) %>% 
  summarise(mean_cpue = mean(value,na.rm=T)) %>% 
  group_by(spp,season) %>% 
  summarise(total_cpue = sum(mean_cpue, na.rm = T))


# Estimate centroid based on top 10 grids with higer abundance
centroid_db <- historic_spp %>% 
  filter(!is.na(label)) %>% 
  group_by(label,index,spp,season,lat,lon) %>% 
  summarise(mean_cpue = mean(value,na.rm=T)) %>% 
  filter(!is.na(mean_cpue)) %>% 
  left_join(total_inter) %>% 
  mutate(
    per_cpue = round((mean_cpue/total_cpue)*100,2)
  ) %>% 
  group_by(spp,label,season) %>% 
  top_n(10,per_cpue) %>% 
  group_by(spp,label,season) %>% 
  summarise(
    lat = mean(lat),
    lon = mean(lon)
  )

# Paper numbers
# centroid_db %>% 
#   select(-lon) %>% 
#   spread(label,lat) %>% 
#   mutate(diff = Today-Reference) %>% 
#   View()


# For alternative color pallet
pal <- wes_palette("Zissou1",100,type = "continuous")

inter_map <- ggplot() +
  geom_tile(data = subset(interpol_result, mean_cpue == 0),
            aes(
              x = lon,
              y = lat
            ),
            fill = "grey90",
            color = "grey90"
  ) + 
  geom_tile(data = subset(interpol_result, mean_cpue > 0),
            aes(
              x = lon,
              y = lat,
              fill = log10(mean_cpue),
              color = log10(mean_cpue)
              # fill = mean_cpue,
              # color = mean_cpue
            )
  ) + 
  geom_point(data = centroid_db,
             aes(
               x = lon,
               y = lat
             ),
             shape = 3,
             color = "white"
  ) + 
  geom_sf(data = land_sf,aes()) +
  labs(x = "Longitude",
       y = "Latitude") +
  # ggsflabel::geom_sf_label_repel(data = land_sf, aes(label= abrev), 
  #                                  size = 4,
  #                                  box.padding = 0.10,
  #                                  hjust = 1) +
  # Viridis option
# scale_fill_viridis("CPUE",
#                    limits = c(0,40),
#                    # limits = c(-5,1), # for log10
#                    direction = -1,
#                     end = 0.9
#                    ) +
# scale_color_viridis("CPUE",
#                     limits = c(0,40),
#                     # limits = c(-5,1), # for log10
#                     direction = -1,
#                     end = 0.9
#                     ) +
# wes anderson option
scale_fill_gradientn("CPUE",
                     colours = pal,
                     # limits = c(0,40) # For cpue
                     limits = c(-5,1) # for log10

) +
  scale_color_gradientn("CPUE",
                        colours = pal,
                        # limits = c(0,40) # for cpue
                        limits = c(-5,1) # for log10
  ) +
  facet_grid(spp~season+label) +
  my_ggtheme_m(map_type = "na", leg_pos = "right",ax_tx_s = 7,ax_tl_s = 10,leg_tx_s = 12) #+
# theme(legend.key.width = unit(3,"line"))


ggsave("interpolation_wa.jpg",
       inter_map,
       path = my_path("R","Figures"),
       width = 9,
       height = 7)

Average proportion

This map shows the aggregated extrapolated value from all three species per State average across the whole study period within each State’s water.

Note: This is intended to be a supplemental figure

data_grid <- historic_spp %>% 
  group_by(state,lat,lon,region,spatial) %>% 
  summarise(sum = sum(value,na.rm= T), .groups = "drop") %>% # sum of all species
  group_by(state,lat,lon,region,spatial) %>% 
  summarise(mean = mean(sum,na.rm= T), .groups = "drop") # average of all years


ggplot(data = subset(data_grid, !is.na(mean))) +
  geom_tile(
    aes(
      x = lon,
      y = lat,
      fill = mean,
      col = mean
    )
  ) +
  geom_sf(data = land_sf, aes()) +
  facet_wrap(~ state  +  region + spatial) +
  labs(x = "Longitude", y = "Latitude") +
  scale_fill_viridis("Mean Interpolation", alpha = 0.8) +
  scale_color_viridis("Mean Interpolation", alpha = 0.8) +
  theme_bw()

Proportion Change

Here we show the average proportion of the interpolation by State and time period. Time periods were arbitrary defined as;

• Early; 1973 to 1984
• Mid; 1985 to 1997
• Late; 1998 to 2011
• Now; 2012 to 2019

Notes: Figure represents the Spring survey. This computation considers the Overlapping of state waters.

Data preparation

total_fited <- historic_spp %>% 
  group_by(year,season,spp,spatial) %>% 
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")


state_fit <- historic_spp %>% 
  group_by(state,year,label,season,spp,spatial) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  # View()
  left_join(total_fited,
            by = c("year","season","spp","spatial")) %>%
  mutate(percentage = state_value/total_value*100) %>% 
  group_by(state,label,season,spp,spatial) %>% 
  summarise(mean_per = round(mean(percentage)),.groups = "drop") %>% 
  #Only show results for spring
  filter(!is.na(label))

Regulatory units map

reg_unit <- state_fit %>% 
  filter(season == "Spring") %>% 
  spread(spatial,mean_per) %>% 
  # Set NA's to 0 because no data found on FP
  mutate(fp = replace_na(fp,0)) %>% 
  # View()
  mutate(
    difference = abs(fp-sw),
    difference = ifelse(difference > 100,100,
                        ifelse(difference < -100,-100,difference)
    )
  ) %>% 
  gather("spatial","mean_per",fp:difference) %>% 
  mutate(spp = gsub(" ","\n",spp),
         spatial = ifelse(spatial == "fp","Fishing ports",
                          ifelse(spatial == "sw","State waters","Difference")
         )
  )


# The plot
map_plot <- land_sf %>% 
  left_join(reg_unit,
            by = "state") %>% 
  filter(spatial != "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Average proportion per State", alpha = 0.8, 
                              breaks = seq(0,35,5),
                              limits=(c(0,35))
                              )+
  facet_grid(spatial~ spp+label) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line")) 


diff_plot <- land_sf %>% 
  left_join(reg_unit,
            by = "state") %>% 
  filter(spatial == "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Difference (%)", alpha = 0.8,option = "C",
                              breaks = seq(0,20,5),
                              limits=(c(0,22)),
                              # labels = c("<-10",seq(-75,75,25),">100")
                              ) +
  facet_grid(spatial~ spp+label) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        strip.background = element_blank(),
        strip.text.x = element_blank()
        ) 


# Cowplot option
p <- ggdraw() +
  # Revenue circular
  draw_plot(map_plot, x = 0, y = 0.2, width = 1, height = 0.8) +
  # Catch circular
  draw_plot(diff_plot, x = 0, y = 0, width = 1, height = 0.4) +
  draw_plot_label(label = c("Latitude", "Longitude"), 
                  size = 18,
                  angle = c(90,0),
                  x = c(0,0.45), 
                  y = c(0.45,0.15)
  )
  

save_plot(
  path = my_path("R","Figures"),
  "proportion_chg_reg.jpg",
  p,
  base_height = 14,
  base_width = 14
)

reg_unit %>% 
  # filter(spatial == "fp") %>%
  select(-fp) %>% 
  spread(label,sw) %>% 
  mutate(
    difference = Today-Reference
    ) %>% 
  gther(Reference:difference)

Seasonal units map

# Get seasonal data
season_fit <- state_fit %>% 
  filter(spatial != "fp") %>%
  spread(season,mean_per) %>%
  mutate(
    Difference = abs(Fall-Spring),
    Difference = ifelse(Difference > 100,100,
                        ifelse(Difference < -100,-100,Difference)
    )
  ) %>% 
  gather("season","mean_per",Fall:Difference) %>% 
  mutate(spp = gsub(" ","\n",spp))

# The plot
map_plot <- land_sf %>% 
  left_join(season_fit,
            by = "state") %>% 
  filter(season != "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Average proportion per State", alpha = 0.8, 
                              breaks = seq(0,15,5),
                              limits=(c(0,15))
                              )+
  facet_grid(season~ spp+label) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line")) 


diff_plot <- land_sf %>% 
  left_join(season_fit,
            by = "state") %>% 
  filter(season == "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Difference (%)", alpha = 0.8,option = "C",
                              breaks = seq(0,5,1),
                              limits=(c(0,5)),
                              # labels = c("<-10",seq(-75,75,25),">100")
                              ) +
  facet_grid(season~ spp+label) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        strip.background = element_blank(),
        strip.text.x = element_blank()
        ) 


# Cowplot option
p <- ggdraw() +
  # Revenue circular
  draw_plot(map_plot, x = 0, y = 0.2, width = 1, height = 0.8) +
  # Catch circular
  draw_plot(diff_plot, x = 0, y = 0, width = 1, height = 0.4) +
  draw_plot_label(label = c("Latitude", "Longitude"), 
                  size = 18,
                  angle = c(90,0),
                  x = c(0,0.45), 
                  y = c(0.45,0.15)
  )
  

save_plot(
  path = my_path("R","Figures"),
  "proportion_chg_seas_.jpg",
  p,
  base_height = 14,
  base_width = 14
)

Area plot

This graph shows the proportion of the interpolation value each State has over time.

Note: This plot is interactive. For ease comparison between States you can;

• One click on any State to remove it from the plot
• Two clicks on any State to isolate it from the plot (other states can then be added by just clicking on them).
• The bottom panel allows you to reduce the time frame

Dynamic aggregated plot

total_fited <- historic_spp %>% 
  group_by(year,region) %>% 
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")

# group by state
state_fit <- historic_spp %>% 
  group_by(state,year,region) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region")) %>%
  mutate(percentage = state_value/total_value*100)

# Plot me
p <- ggplot(state_fit) +
  geom_area(
    aes(
      x = year,
      y = round(percentage),
      fill = state
    )
  ) +
  ylab("Percentage (%)") +
  # viridis::scale_fill_viridis(discrete = T, alpha = 0.8) +
  scale_fill_manual(values = state_pallet) +
  MyFunctions::my_ggtheme_p() +
  facet_wrap(~region, ncol = 1) +
  theme_dark()

ggplotly(p,
         dynamicTicks = TRUE) %>% 
  layout(hovermode = "x") %>%
  rangeslider()

Area plot (running mean)

This graph shows the proportion of each State smoothed over a *10 years running mean**. It allows to better see increasing and decreasing trends.

Dynamic aggregated plot

Note: This plot is also interactive and thus, follows the same options of the previous one.

# group by state
state_fit <- historic_spp %>% 
  group_by(state,year,region,.groups = "drop") %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region")) %>%
  mutate(percentage = state_value/total_value*100) %>% 
  group_by(state,region) %>% 
  mutate(RMean = zoo::rollmean(x = percentage, 
                            10,
                            fill = NA,
                            na.rm=T)
    ) %>%  
  left_join(state_lat)

# Plot me
p <- ggplot(state_fit) +
  geom_area(
    aes(
      x = year,
      y = round(RMean),
      fill = state
      # fill = order
    )
  ) +
  ylab("10 yrs running average (%)") +
  scale_fill_manual(
    "State",
    values = state_pallet,
    # labels = state_fit %>% arrange(order) %>%  pull(state) %>% unique()
  ) +
  MyFunctions::my_ggtheme_p() +
  facet_wrap(~region, ncol = 1) +
  theme_dark()

suppressWarnings(
ggplotly(p,
         dynamicTicks = TRUE) %>% 
  layout(hovermode = "x") %>% 
  # add_trace() %>% 
  rangeslider()
)

Aggregated Proportional change

Paper figure: Per Species Proportional change

Old code

Fishing ports WPI

We used Global Fishing Watch’s database on Anchorages, Ports and Voyages.

The Global Fishing Watch anchorages dataset is a global database of anchorage locations where vessels congregate. The dataset contains over 160,000 individual anchorage locations, which are further grouped into nearly 32,000 ports when applicable.

# Get Global Fishig Watch data
gfw_data <- my_path("D","Spatial",name = "gfw_portdata.csv", read = T) %>% 
  filter(iso3 == "USA") 

# Visualize data (o.k!)
# ggplot(gfw_data) +
#   geom_point(
#     aes(
#     x = lon,
#     y = lat
#   )
#   )

# Convert data to sf for mergig
gfw_sf <- st_as_sf(gfw_data,
             coords = c("lon", "lat"),
             crs = 4326)

# Merge data with state waters to,locate ports 
port_sf <- st_join(gfw_sf,state_sf) %>% 
  filter(!is.na(jurisdicti)) # remove points outside stat waters

ggplot(port_sf)  +
  geom_sf(aes(color = distance_from_shore_m)) +
  theme_dark()

So…Looks like some ports are “far from shore” which I don’t really know what it means. For now, I am only selecting those ports that ar 0 meter from shore.

There are still too many ports in each state (?). Specially New Jersey! One approach is to use NOAA’s landing-per-port data to filter out…

4.2 World Port Index

In this approach we only use the ports categorized within the WPI and with distance_from_shore_m == 0 and at_dock == TRUE.

• at_dock == Anchorages within a 450 meter buffer around a combined land shape file are considered at dock.

Historic quota Vs distributions

This section of the results explores the differences between the historic distribution of te stock and the historic quota allocation. We collected the proportion of the total quota that each State gets for each species.

Right now the analysis only covers black sea bass, summer flounder and scup. But we can go speceis by species to see which ones have their quota allocated by state to include them in the analysis.

• Black Sea Bass; Centropristis striata

Note: We are using the new allocations based on the stock’s most recent biomass distribution*

• Summer flounder

• Scup Summer period

For State-level managed species see the Atlantic States Marine Fisheries Comission website and go species-by-species.

quota_allocation <- state_lat %>% 
  mutate(
    "Centropristis striata" = c(
      0.0040, #ME
      0.0040, #NH
      0.1564, #MA
      0.1323, #RI
      0.0367, #CT
      0.0857, #NY
      0.2010, #NJ
      0.0411, #DE
      0.0888, #MD
      0.1614, #VA
      0.0888 #NC
    ),
    "Paralichthys dentatus" = c(
      0.0004756, #ME 
      0.0000046, #NH 
      0.0682046, #MA 
      0.1568298, #RI 
      0.0225708, #CT 
      0.0764699, #NY 
      0.1672499, #NJ 
      0.0001779, #DE 
      0.0203910, #MD
      0.2131676, #VA
      0.2744584 #NC
    ),
    "Stenotomus chrysops" = c(
      0.0012101, #ME
      0.0000000, #HN
      0.2158729, #MA
      0.5619456, #RI
      0.0315399, #CT
      0.1582466, #NY
      0.0291667, #NJ
      0.0000000, #DE
      0.0001190, #MD
      0.0016502, #VA
      0.0002490 #NC
    )
  ) %>% 
  gather("spp","quota",4:6)

# Double check they all add to 1...
quota_allocation %>% 
  group_by(spp) %>% 
  summarise_at(vars(quota),sum)


quota_allocation %>% 
  arrange(order) %>% 
  select(-order) %>% 
  spread(spp,quota) %>% 
  kable()

Get data ready

total_fited <- historic_spp %>% 
  group_by(year,region,spp) %>% 
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")

# group by state
state_fit <- historic_spp %>% 
  group_by(state,year,region,spp) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region","spp")) %>%
  mutate(percentage = state_value/total_value*100) %>% 
  #Only show results for spring
  filter(str_detect(region,"Spring"))

quota_df <- quota_allocation %>% 
  left_join(historic_spp) %>% 
  group_by(state,abrev,year,region,spp,quota) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region","spp")) %>%
  mutate(Distribution = state_value/total_value*100,
         Historic = quota*100) %>% 
  group_by(state,quota,spp) %>% 
  mutate(RMean = zoo::rollmean(x = Distribution, 
                            10,
                            fill = NA,
                            na.rm=T)
    ) %>% 
  #Only show results for spring
  filter(str_detect(region,"Spring"),
         !is.na(RMean) # Remove NAs from rollmean
         ) %>%  
  left_join(state_lat) %>% 
  filter(
    # spp =="Centropristis striata"
    spp %in% quota_allocation$spp
  )

head(quota_df) %>% 
  kable()

Exploring different ways to show the results

Difference line plot

This version shows the difference between the Historic quota allocation and the distribution proportion of the stock:

• diff < 0 = the distribution proportion is larger than the allocated quota
• diff > 0 = the distribution proportion is smaller than the allocated quota

# Option one

quota_df %>%
  mutate(diff = Historic- Distribution,
         diff_rmean = RMean- Distribution) %>% 
  gather("type","value",diff:diff_rmean) %>% 
  # View()
  ggplot() +
  geom_line(
    aes(
      x = year,
      y = value,
      color = order
    )
  ) +
  scale_color_manual("State",
    values = state_pallet,
    labels = quota_df %>% arrange(order) %>%  pull(state) %>% unique()
    ) +
  facet_grid(type~spp, 
             scales = "free") +
  theme_dark()

Linetype plot

Here, we plot the distributional quota (solid lines) with each State’s allocation with dashed (—-) lines

# Option one

quota_df %>% 
  gather("level","quota",Distribution:RMean) %>% 
  # mutate(diff = quota_per- percentage) %>% 
  # View()
  ggplot() +
  geom_line(
    aes(
      x = year,
      y = quota,
      color = order
    )
  ) +
  labs(x = "Year",y = "Quota (%)") +
  scale_fill_manual("State",
    values = state_pallet,
    labels = quota_df %>% arrange(order) %>%  pull(state) %>% unique()
    ) +
  scale_color_manual("State",
    values = state_pallet,
    labels = quota_df %>% arrange(order) %>%  pull(state) %>% unique()
    ) +
  scale_linetype("Quota") +
  facet_grid(spp ~ level)+
  theme_dark()

Efficiency index

The idea here is create and efficiency index (Danger!) that computes the alignment between the distribution and the actual allocation. The index is simply;

\[ei = \frac{HistoricQuota}{DistributionProportion}\] So, in cases where ei = 1, the quota is aligned with the distribution; when ei > 1 then the historic quota is more than the stock’s State’s distribution, contrarily, ei < 1 means that the historic quota is less than what the State currently has.

Note: - There are some crazzy outliers that have been removed, for now… - Bottom plot representing index as a Rmean

As a normal year to year

# Option one

eff_index <- quota_df %>% 
  # filter(
  #   spp =="Centropristis striata"
  # ) %>% 
  mutate(index = Distribution/Historic,
         index_rmean = RMean/Historic) %>%
  gather("level","index",index:index_rmean) %>%
  filter(index < 5) # there are some craaaazy outliers
  # View()
  
ggplot(eff_index) +
  geom_line(
    aes(
      x = year,
      y = index,
      color = state
    ),
  ) +
  labs(x = "Year",y = "Efficiency index") +
  scale_fill_manual("State",
    values = state_pallet,
    labels = quota_df %>% arrange(order) %>%  pull(state) %>% unique()
    ) +
  scale_color_manual("State",
    values = state_pallet,
    labels = quota_df %>% arrange(order) %>%  pull(state) %>% unique()
    ) +
  facet_grid(level~spp,
             scales = "free")+
  theme_dark()

NOAAs port approach

Dynamic figure: Aggregated Proportional change per survey

total_fited <- historic_spp %>% 
  group_by(year,region,spatial) %>% 
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")


state_fit <- historic_spp %>% 
  group_by(state,year,region,spatial) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region","spatial")) %>%
  mutate(percentage = state_value/total_value*100) %>% 
  left_join(periods,
            by = "year") %>% 
  group_by(state,order,label,region,spatial) %>% 
  summarise(mean_per = round(mean(percentage)),.groups = "drop")

# check percentages
# state_fit %>% 
#   group_by(period) %>% 
#   summarise(sum(mean_per))


p <- land_sf %>% 
  left_join(state_fit,
            by = "state") %>% 
  ggplot() +
  # geom_sf(aes(fill = mean_per)) +
  geom_sf(aes(fill = mean_per, text = paste(state, mean_per,"% of proportion"))) +
  viridis::scale_fill_viridis("Average proportion per State", alpha = 0.8) +
  facet_wrap(~ region + label +spatial, 
             nrow = 2) +
  labs(x = "Latitude", 
       y = "Longitude") +
  my_ggtheme_p(facet_tx_s = 12,
               leg_pos = "bottom") +
  theme(legend.key.width = unit(1,"line")) +
  theme_dark()

ggplotly(p,
         tooltip = "text",
         dynamicTicks = FALSE) %>% 
    style(hoverlabel = list(bgcolor = "white"), hoveron = "fill")

total_fited <- historic_spp %>%
  group_by(year,region,spp,spatial) %>%
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")


state_change <-historic_spp %>% 
  group_by(state,year,region,spp,spatial) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region","spp","spatial")) %>%
  mutate(percentage = state_value/total_value*100,
         label = ifelse(year >= 1980 & year <= 2001,"Reference",
                         ifelse(year > 2001,"Today",NA)
                         )
         ) %>% 
  group_by(state,label,region,spp,spatial) %>% 
  summarise(mean_per = round(mean(percentage)),.groups = "drop") %>% 
  # Only show results for spring
  filter(str_detect(region,"Spring"),
         !is.na(label)
         ) %>% 
  select(-region) %>% 
  spread(label,mean_per) %>% 
  mutate(change = (Today-Reference)/(Today)*100,
         change = ifelse(change > 100,100,
                         ifelse(change < -100,-100,change)
         )
  ) %>%
  select(-4,-5) %>% 
  spread(spatial,change) %>% 
  mutate(
    #difference = (fp-sw)/((fp+sw)/2)*100,
    difference = abs(fp-sw)#,
         # difference = ifelse(difference > 100,100,
         #                     ifelse(difference < -100,-100,difference)
         # )
  ) %>% 
  gather("spatial","mean_per",fp:difference) %>% 
  mutate(
         spatial = ifelse(spatial == "fp","Fishing ports",
                          ifelse(spatial == "sw","State Waters","Difference")
         ),
         mean_per= replace_na(mean_per,0)# for those 0/0
  )
   

# The plot
prop_chng_map <- land_sf %>% 
  left_join(state_change,
            by = "state") %>% 
  filter(spatial != "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Change in proportion today\nrelative to 1980-2001\n(%)", alpha = 0.8,
                              limits = c(-100,100),
                              breaks = seq(-100,100,25),
                              labels = c("<-100",seq(-75,75,25),">100")
                              ) +
  facet_grid(spatial~ spp) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        legend.title.align=0.5) 

# ggsave(filename = "temp_proportion_chng.jpg",
#          plot = prop_chng_map,
#          path = my_path("R","Figures"),
#          width = 10,
#          height = 10
#   )


#### Option with differences

prop_chng_map_diff <- land_sf %>% 
  left_join(state_change,
            by = "state") %>% 
  filter(spatial == "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Change in proportion today\nrelative to 1980-2001\n(%)", alpha = 0.8,
                              option = "C"
                              limits = c(-100,100),
                              breaks = seq(-100,100,25),
                              labels = c("<-100",seq(-75,75,25),">100")
                              ) +
  facet_grid(spatial~ spp) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        legend.title.align=0.5,
        strip.background = element_blank(),
        strip.text.x = element_blank())


# Cowplot option
p <- ggdraw() +
  # Revenue circular
  draw_plot(prop_chng_map, x = 0, y = 0.25, width = 1, height = 0.65) +
  # Catch circular
  draw_plot(prop_chng_map_diff, x = 0.101, y = 0, width = 0.798, height = 0.4) +
  draw_plot_label(label = c("Latitude", "Longitude"), 
                  size = 16,
                  angle = c(90,0),
                  x = c(0.1,0.45), 
                  y = c(0.45,0.11)
  )
    

save_plot(
  path = my_path("R","Figures"),
  "proportion_chg_diffN.jpg",
  p,
  base_height = 14,
  base_width = 14
)

Proportion change realtive to past

total_fited <- historic_spp %>%
  group_by(year,season,spp,spatial) %>%
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")


state_change <-historic_spp %>% 
  group_by(state,year,season,spp,spatial) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region","spp","spatial")) %>%
  mutate(percentage = state_value/total_value*100,
         label = ifelse(year >= 1980 & year <= 2001,"Reference",
                         ifelse(year > 2001,"Today",NA)
                         )
         ) %>% 
  group_by(state,label,region,spp,spatial) %>% 
  summarise(mean_per = round(mean(percentage)),.groups = "drop") %>% 
  # Only show results for spring
  filter(str_detect(region,"Spring"),
         !is.na(label)
         ) %>% 
  select(-region) %>% 
  spread(label,mean_per) %>% 
  mutate(change = (Today-Reference)/(Today)*100,
         change = ifelse(change > 100,100,
                         ifelse(change < -100,-100,change)
         )
  ) %>%
  select(-4,-5) %>% 
  spread(spatial,change) %>% 
  mutate(
    #difference = (fp-sw)/((fp+sw)/2)*100,
    difference = abs(fp-sw)#,
         # difference = ifelse(difference > 100,100,
         #                     ifelse(difference < -100,-100,difference)
         # )
  ) %>% 
  gather("spatial","mean_per",fp:difference) %>% 
  mutate(
         spatial = ifelse(spatial == "fp","Fishing ports",
                          ifelse(spatial == "sw","State Waters","Difference")
         ),
         mean_per= replace_na(mean_per,0)# for those 0/0
  )
   

# The plot
prop_chng_map <- land_sf %>% 
  left_join(state_change,
            by = "state") %>% 
  filter(spatial != "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Change in proportion today\nrelative to 1980-2001\n(%)", alpha = 0.8,
                              limits = c(-100,100),
                              breaks = seq(-100,100,25),
                              labels = c("<-100",seq(-75,75,25),">100")
                              ) +
  facet_grid(spatial~ spp) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        legend.title.align=0.5) 

# ggsave(filename = "temp_proportion_chng.jpg",
#          plot = prop_chng_map,
#          path = my_path("R","Figures"),
#          width = 10,
#          height = 10
#   )


#### Option with differences

prop_chng_map_diff <- land_sf %>% 
  left_join(state_change,
            by = "state") %>% 
  filter(spatial == "Difference") %>% 
  ggplot() +
  geom_sf(aes(fill = mean_per)) +
  viridis::scale_fill_viridis("Change in proportion today\nrelative to 1980-2001\n(%)", alpha = 0.8,
                              option = "C"
                              limits = c(-100,100),
                              breaks = seq(-100,100,25),
                              labels = c("<-100",seq(-75,75,25),">100")
                              ) +
  facet_grid(spatial~ spp) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        legend.title.align=0.5,
        strip.background = element_blank(),
        strip.text.x = element_blank())


# Cowplot option
p <- ggdraw() +
  # Revenue circular
  draw_plot(prop_chng_map, x = 0, y = 0.25, width = 1, height = 0.65) +
  # Catch circular
  draw_plot(prop_chng_map_diff, x = 0.101, y = 0, width = 0.798, height = 0.4) +
  draw_plot_label(label = c("Latitude", "Longitude"), 
                  size = 16,
                  angle = c(90,0),
                  x = c(0.1,0.45), 
                  y = c(0.45,0.11)
  )
    

save_plot(
  path = my_path("R","Figures"),
  "proportion_chg_diffN.jpg",
  p,
  base_height = 14,
  base_width = 14
)

Aggregated Proportional change

total_fited <- historic_spp %>%
  group_by(year,region,spatial) %>%
  summarise(total_value = sum(value,na.rm=T),.groups = "drop")


state_change <- historic_spp %>% 
  group_by(state,year,region,spatial) %>% 
  summarise(state_value = sum(value,na.rm= T), .groups = "drop") %>% 
  left_join(total_fited,
            by = c("year","region","spatial")) %>%
  mutate(percentage = state_value/total_value*100,
         label = ifelse(year >= 1980 & year <= 2001,"reference",
                         ifelse(year > 2001,"today",NA)
                         )
         ) %>% 
  # left_join(periods,
            # by = "year") %>% 
  group_by(state,label,region,spatial) %>% 
  summarise(mean_per = round(mean(percentage)),.groups = "drop") %>% 
  #Only show results for spring
  filter(#str_detect(region,"Spring"),
         !is.na(label)
         ) %>% 
  # select(-region) %>% 
  # Estimate percentage change
  spread(label,mean_per) %>% 
  mutate(raw_change = (today-reference)/(today)*100,
         change = ifelse(raw_change > 100,100,
                         ifelse(raw_change < -100,-100,raw_change)
         )
  ) %>%
  select(-reference,-today) 


#### Average
state_change %>% 
  mutate(
    wl = ifelse(change < 0, "L",
                ifelse(change == 0, "NC","W"))
  ) %>% 
  # View()
  # group_by(wl) %>% 
  group_by(region,spatial,wl) %>% 
  summarise(
    states = paste(unique(state),collapse = ", "),
    mean(raw_change),
    sd(raw_change)
  ) %>% 
  write_csv("../Results/state_relative_chng.csv")
  View()



  # Estimate spatial differences in percentage change
spatial_diff <- state_change %>% 
  spread(spatial,change) %>% 
  mutate(spatial_diff = abs(fp-sw)) %>% 
  select(state, category = region, diff = spatial_diff)

# Estimate seasonal differences in percentage change
season_diff <- state_change %>% 
  spread(region,change) %>%  
  mutate(season_diff = abs(`Northeast US Fall`-`Northeast US Spring`)) %>% 
  select(state,category = spatial,diff = season_diff)

differences <- season_diff %>% 
  bind_rows(spatial_diff) %>% 
  mutate(
    label = ifelse(category == "fp","Fishing ports",
                          ifelse(category == "sw","State Waters",category)
         )
  )

# For in text differences
differences %>% 
  group_by(category) %>%
  summarise_at(vars(diff),
    c(mean = mean,
      sd =sd,
      min = min,
      max= max))


# The main plot
prop_chng_map <- land_sf %>% 
  left_join(state_change,
            by = "state") %>% 
  # filter(spatial != "Difference" & cat != "difference_reg") %>% 
  ggplot() +
  geom_sf(aes(fill = change)) +
  viridis::scale_fill_viridis("Change in proportion today\n relative to 1973-1984\n(%)", alpha = 0.8,
                              limits = c(-100,100),
                              breaks = seq(-100,100,25),
                              labels = c(seq(-100,75,25),">100")
                              ) +
  facet_grid(spatial~ region) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        legend.title.align=0.5,
        strip.text.y = element_blank(),
        strip.background = element_blank(),
        ) 

ggsave(filename = "proportion_chng_agg.png",
         plot = prop_chng_map,
         path = my_path("R","Figures"),
         width = 10,
         height = 10
  )


#### Option with differences

# Season differencce

spat_diff <- land_sf %>% 
  left_join(differences,
            by = "state") %>% 
  filter(!category %in% c("fp","sw")) %>% 
  ggplot() +
  geom_sf(aes(fill = diff)) +
  viridis::scale_fill_viridis("Differences", 
                              alpha = 0.8,
                              option = "C",
                              limits = c(0,50),
                              breaks = seq(0,50,10)
                              ) +
  facet_wrap(~ label,
             ncol = 2) +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "bottom",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(legend.key.width = unit(4,"line"),
        legend.title.align=0.5,
        strip.background = element_blank(),
        strip.text.x = element_blank())

ggsave(filename = "spat_diff.png",
         plot = spat_diff,
         path = my_path("R","Figures"),
         width = 8,
         height = 6
  )

# seasonal
seas_diff <- land_sf %>% 
  left_join(differences,
            by = "state") %>% 
  filter(category %in% c("fp","sw")) %>% 
  ggplot() +
  geom_sf(aes(fill = diff)) +
  viridis::scale_fill_viridis("Differences", 
                              alpha = 0.8,
                              option = "C",
                              limits = c(0,50),
                              breaks = seq(0,50,10)
                              ) +
  facet_wrap(~ label,
             ncol = 1,
             strip.position="right") +
  labs(x = "", 
       y = "") +
  my_ggtheme_p(facet_tx_s = 20,
               leg_pos = "",
               axx_tx_ang = 45,
               ax_tx_s = 12,
               ax_tl_s = 18,
               hjust = 1) +
  theme(strip.background = element_blank(),
        # strip.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.line.y = element_blank()
        )

ggsave(filename = "seas_diff.png",
         plot = seas_diff,
         path = my_path("R","Figures"),
         width = 8,
         height = 6
  )

### SRAH LANDING DATA
# read_excel_allsheets <- function(filename, tibble = TRUE) {
#     # I prefer straight data.frames
#     # but if you like tidyverse tibbles (the default with read_excel)
#     # then just pass tibble = TRUE
#     sheets <- readxl::excel_sheets(filename)[2:8]
#     x <- lapply(sheets, function(X) readxl::read_excel(filename, sheet = X))
#     if(!tibble) x <- lapply(x, as.data.frame)
#     names(x) <- sheets
#     x
# }
# Get data path
# ports_path <- my_path("D","NOAA", name = "SSmith_2015-2021 Landings_JUN 2021.xlsx")

# Load ports
# ports_data <- bind_rows(read_excel_allsheets(ports_path)) %>%
#   clean_names() %>%
#   mutate_if(is.character,
#               ~ str_to_lower(.)) %>%
#     filter(species_name %in% c("summer flounder","black sea bass","scup"))